package scu.maqiang.homogenization;

import java.util.Arrays;

import scu.maqiang.fes.*;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.*;

public class PossionCell3DT4 {

	public static void main(String[] args) {
		long startTime = System.nanoTime();
        int N = 1, label = 10;
        double ep = 1.0 / N;
		Mesh3T4 mesh = new Mesh3T4().cube(40, 40, 40);
        ScalarFunc cellRegionFunc = (xyz, llable, param) -> {
            return (Math.abs(xyz[0] - 0.5) < 0.25 && Math.abs(xyz[1] - 0.5) < 0.25 && Math.abs(xyz[2] - 0.5) < 0.25)? 1.0:-1.0;
        };
        
        mesh.setDomainLabel(cellRegionFunc, null, label);
//        mesh.displayElements();
        
        FES3T41 fs = new FES3T41(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        ScalarFunc regionCoef = (xyz, mylabel, param) -> mylabel == label ? 0.01 : 1.0;
        ScalarFunc[] gcf = {regionCoef};
        fs.assembleStiff(gcf, null, BVPType.COMMON, A);
        
//        A.sort();
//        System.out.println(A);
		System.out.println("Size Of A = " + A.getNNZ());
		double[][] RHS = new double[3][fs.GetNdof()];
		fs.assembleSourceDerivative(gcf, null, BVPType.COMMON, RHS);
		fs.applyBC_MBN(A, Direct.All, 1, 2, 3, 4, 5, 6);
		fs.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 3, 4, 5, 6);
		fs.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 3, 4, 5, 6);
		fs.applyBC_MBN(RHS[2], Direct.All, 0.0, 1, 2, 3, 4, 5, 6);
		
		double[][] x = new double[3][fs.GetNdof()];
		NewIterSSolver newSolver = new NewIterSSolver(A);
		newSolver.PCGSSOR(RHS, x, 1.5e0, 1);
//		mesh.toTecplot("CellSolution.dat", x);
		
		Arrays.fill(x[0], 0.0);
		newSolver.PCGSSOR(RHS[0], x[0], 1.5e0, 1);
		Arrays.fill(x[1], 0.0);
		newSolver.PCGSSOR(RHS[1], x[1], 1.5e0, 1);
		Arrays.fill(x[2], 0.0);
		newSolver.PCGSSOR(RHS[2], x[2], 1.5e0, 1);
		
		mesh.toTecplot("CellSolution3D.dat", x);
//		double m = fespace.integrate(coef, QuadType.Gauss1);
//		System.out.println("Volume Average = " + m);
//		double[][] m1 = fespace.computeHomogeniztionHeat(coef, N1, N2);		
//		System.out.println("Homogenization Value = \n" + MMatrix.toString(m1));
		
		long endTime = System.nanoTime();
		System.out.println("Elapse Time: " + (endTime - startTime) / 1e9 + " s");
	}
}